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Chapter  1:  Introduction 


Since  the  first  surveys  were  conducted  in  1982,  numerous  GPS  applications 
have  been  implemented  due  to  its  high  degree  of  accuracy.  Whether  military  or 
civilian  in  nature,  GPS  has  redefined  navigation  on  Earth  as  well  as  in  space.  It  is 
being  used  in  orbit  determination  and  is  beginning  to  play  an  important  role  in 
spacecraft  attitude  determination.  One  of  the  limitations  to  GPS,  however,  is  its 
dependence  on  line-of-sight  signal  propagation  between  the  satellite  constellation  and 
the  receiver.  In  certain  terrestrial  applications  where  an  imobstructed  view  of  the  sky 
is  not  available,  pseudolites  have  been  implemented  to  fill  the  void. 

The  National  Aeronautics  and  Space  Administration  -  Johnson  Space  Center 
(NAS  A- JSC),  in  conjunction  with  the  University  of  Texas  at  Austin  (UT)  Center  for 
Space  Research  (CSR),  is  researching  the  applicability  of  pseudolites  (GPS-like 
signal  propagators)  to  create  a  body  fixed  navigation  reference  fi-ame  for  elose 
proximity  space  vehicle  operations  [1].  More  specifically,  the  International  Space 
Station  (ISS)  could  benefit  fi'om  an  autonomous  navigation  system  based  on  GPS 
signal  theory  because  docking  requirements,  close  proximity  formation  flight,  and 
Crew  Return  Vehicle  (CRV)  flight  paths  all  require  accurate  navigation.  However, 
these  maneuvers  could  easily  lose  the  required  line-of-sight  to  the  GPS  eonstellation 
when  performed  elose  to  the  large  space  station.  Before  a  GPS-independent  system 
can  be  implemented,  however,  it  is  important  to  understand  the  differences  between 
GPS  and  pseudolite  signals. 
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This  paper  first  explores  the  navigation  solution  required  to  pinpoint  the 
location  of  a  static  pseudolite  receiver  in  an  indoor  JSC  Navigation  Systems  and 
Technology  Lab  (NSTL)  using  rmbiased  pseudoranges  formed  fi-om  code  phase 
signals.  It  then  applies  a  carrier  phase  smoothing  technique  to  form  a  more  accurate 
post-processed  data  solution  because  the  pseudolite  signals  are  propagated  extremely 
close  to  the  receiver  in  comparison  to  most  GPS  applications,  and  the  receiver  is 
subjected  to  the  noisy  conditions  inherent  in  the  NSTL.  The  close  proximity  situation 
presented  in  the  lab  is  sirhilar  to  what  will  be  expected  in  the  ISS  reference  fi'ame  in 
terms  of  signal  to  noise  ratios  and  the  unavailability  of  GPS. 

The  initial  goal  of  this  work  was  to  achieve  a  static  positioning  accuracy  of 
five  meters  (a  value  slightly  larger  than  predicted  code  phase  accuracy  due  to  lab 
noise)  in  simulation.  Once  this  was  achieved,  the  same  algorithms  were  to  be  applied 
to  pseudolite  data  taken  in  the  NSTL.  The  results  of  this  research  will  be  used  as  a 
baseline  for  future  relative  and  kinematic  position  fixes  that  can  in  turn  be  used  to 
derive  relative  navigation  and  collision  avoidance  algorithms  in  the  ISS  reference 
fi'ame. 

Incorporating  GPS  and  pseudolite  reception  capability  into  the  same  receiver 
could  eventually  provide  a  cost  effective  means  to  switch  fi’om  a  GPS  navigation 
solution  to  an  autonomous  pseudolite  navigation  system  when  approaching  ISS. 
Furthermore,  such  a  receiver  could  help  automate  close  proximity  operations  and 
reduce  pilot  workload.  Altogether,  a  combined  GPS  and  pseudolite  navigation 

system  would  be  a  valuable  means  to  separate  future  ISS  traffic. 
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Chapter  2:  Pseudolite,  Hardware,  and  Facility  Description 


Pseudolites  are  currently  used  as  a  means  to  augment  the  GPS  satellite 
constellation.  They  are  mainly  used  when  it  is  desirable  to  obtain  GPS  signals  in 
areas  where  the  required  lines  of  sight  to  the  sateUites  are  not  readily  available  or 
where  an  increased  number  of  signals  is  required.  PseudoUtes  are  devices  that  can 
transmit  GPS-like  signals,  and  many  current  applications  include  transceiving  GPS 
signals  at  a  fixed  location  on  the  Earth's  surface  into  areas  'with  a  hindered  view  of  the 
sky.  Such  appUcations  include  deep  pit  mining,  precision  farming,  and  in  future 
aircraft  precision  landing  systems  known  as  local  and  wide  area  augmentation  [2]. 


Figure  1:  Pseudolite  Transceiver  Apparatus  in  NSTL 

All  of  the  pseudolite  data  required  for  this  project  was  obtained  in  the  NSTL, 

or  high  bay  lab  at  JSC.  The  lab’s  original  configuration  had  six.  IntegriNautics 
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IN200C  Pseudolites  (Figure  1)  positioned  so  that  one  was  in  each  of  the  four  top 
comers  of  the  lab  and  two  were  overhead  (Figure  2).  Later  tests  used  a  four- 
pseudolite  configuration  because  two  pseudohtes  were  removed  for  closer  study  at 
CSR.  Each  pseudohte  points  toward  a  reference  marker  located  on  a  table  near  the 
center  of  the  bay.  The  marker  represents  the  (0, 0, 0)  origin  in  a  lab  fixed  coordinate 
system  defined  by  a  theodolite  survey.  The  pseudolites  are  not  moved  during  or  in- 
between  tests  so  that  their  positions  in  the  NSTL  remain  constant. 


Figure  2:  NSTL  High  Bay  Layout  (Not  to  Scale) 

The  pseudolites  used  in  this  research  are  somewhat  simple  compared  to  the 
complexities  involved  with  actual  GPS  satellites.  Since  their  positions  are  always 
fixed,  they  do  not  send  a  navigation  data  message.  Therefore,  only  a  priori 
knowledge  of  their  locations  in  the  lab  coordinate  frame  is  required  for  data 
processing.  They  only  transmit  the  LI  (1575.42  MHz,  19.0  cm)  GPS-like  signal 
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containing  code  and  carrier  phase  data.  Accurate  pseudorange  data,  however,  must 
be  manipulated  from  a  code  phase,  as  will  be  discussed  in  Chapter  3. 

The  IN200C  pseudolite  has  programmable  power  attenuation  and  an  external 
amplification  up  to  +26  dBm  (1  Watt)  [3].  This  is  an  important  feature  because  of  the 
high  levels  of  noise  inherent  in  the  lab  environment.  It  is  important  to  keep  a  high 
signal  to  noise  ratio  (SNR)  while  being  careful  not  to  cross  correlate  the  signals  as 
they  enter  the  receiver  [4].  All  six  pseudolites  in  the  lab  are  referenced  to  a  single 
oseillator,  and  the  line  biases  between  the  oscillator  and  each  pseudolite  have  been 
minimized  since  each  of  the  six  lines  are  approximately  the  same  length  [5].  While 
each  pseudolite  theoretically  ticks  at  the  same  transmit  time,  the  true  time  is  unknown 
to  the  receiver  because  a  navigation  message  containing  a  "GPS  time"  is  not  sent. 

The  reeeiver  used  in  this  research  was  the  Mitel  Architect  "Ship  Channel" 

receiver.  This  GPS  receiver  has  programmable  capability  through  software  that  was 

modified  to  track  the  NSTL  pseudolites.  Data  is  retrieved  in  the  Receiver 

Independent  Exchange  Format  (RINEX)  [6].  The  software  is  compiled  and  uploaded 

to  the  receiver,  and  then  RINEX  data  is  downloaded  and  saved  to  an  external 

computer.  In  this  application,  the  data  corresponds  to  the  static  antenna  position  over 

the  NSTL  lab  origin  on  the  workbench  (Figure  3).  These  RINEX  files  contain 

pseudorange,  carrier  phase,  Doppler,  and  code  phase  data.  The  code  phase  data  has  a 

1 .023e6  chipping  rate  that  must  be  manipulated  into  a  pseudorange  since  the 

transmitted  pseudoranges  are  inaccurate  without  correct  timing  data.  The  Mitel 

Architect  can  accept  a  navigation  data  packet,  so  it  is  possible  for  this  receiver  to 
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receive  real  GPS  data  if  required.  However,  new  software  must  be  loaded  into  the 
receiver  to  switch  between  GPS  mode  and  pseudolite  mode. 


Figure  3:  Static  Workbench  Layout 

The  tables,  equipment,  walls  and  large  metal  high  bay  door  in  the  NSTL 
produce  many  opportunities  for  multipath  to  occur  while  taking  data.  Because  the 
pseudolites  are  so  close  to  the  receiver  (approximately  7  m),  there  is  also  the 
possibility  of  the  signal  not  sufficiently  spreading  before  reaching  the  receiver 
antenna,  resulting  in  delayed  lock  or  loss  of  signals  [7]  [8].  While  difficult  to 
measure,  lab  multipath  requires  careful  consideration  during  data  interpretation. 
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Chapter  3:  Pseudoranges  from  Code  Phases 


Most  of  the  initial  work  for  this  project  required  a  basic  understanding  of  the 
NSTL  setup  and  operation  as  well  as  an  understanding  of  how  the  Mitel  Architect 
receiver  handles  pseudolite  data.  Using  data  taken  previously  at  the  NSTL,  a  Matlab 
script  was  written  to  read  the  RINEX  observables  produced  from  the  Architect.  From 
this  script,  raw  data  plots  could  be  fisrmed  showing  pseudorange,  carrier  phase, 
Doppler,  and  code  phase  data  versus  time.  The  following  discussion  illustrates  how 
the  raw  pseudoranges  produced  by  the  receiver  were  inaccurate  and  how  true 
pseudorange  values  were  calculated  from  the  code  phase  data. 

An  initial  study  of  these  plots  showed  that  a  time  tag  problem  existed  because 
the  receiver  clock  reset  itself  every  14  seconds.  An  inquiry  to  IntegriNautics 
explained  the  problem  [9].  Although  the  pseudolites  do  not  transmit  relevant 
navigation  data,  they  transmit  navigation  packets  in  the  standard  ICD-200  format  so 
receivers  can  still  lock  onto  their  signals.  Because  a  pseudolite  is  not  on-orbit  like 
normal  GPS  satellites  and  because  it  does  not  know  the  correct  GPS  time,  the 
navigation  message  is  useless  in  terms  of  ephemeris  and  timing  data.  The  receiver 
software  was  written  to  reject  an  invalid  GPS  navigation  message,  which  was  causing 
the  software  clock  to  reset.  The  firmware  was  modified  so  the  receiver  would  ignore 
the  contents  of  the  pseudolite  data  message  altogether.  This  eliminated  the  resetting 
time  tag  problem  in  future  data  sets. 
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Pseudolite 

PRN 

X  Coordinate 

Y  Coordinate 

Z  Coordinate 

Distance 

A 

11 

-0.0446039 

-0.7731419 

-6.7365494 

6.781 

B 

8 

-0.0303306 

3.3911386 

-7.0137452 

7.791 

C 

30 

-4.1225911 

-1.9488081 

-3.2345285 

5.591 

D 

20 

5.4105458 

-2.3215789 

-3.3471627 

6.773 

E 

28 

2.9573544 

5.1741947 

-4.1242551 

7.248 

F 

13 

-4.2642563 

4.4336031 

-2.5139543 

6.645 

Average  Distance;  6.805 

Table  1:  Pseudolite  Phase  Center  Locations  (Meters  from  Lab  Origin) 


In  terms  of  the  NSTL  setup,  a  theodolite  survey  was  taken  to  accurately 
determine  the  pseudolite  phase  center  coordinates  in  the  lab  reference  frame  prior  to 
acquiring  any  relevant  data  (Table  1).  The  pseudolites  were  not  moved  again  after 
the  survey.  Also  prior  to  taking  data,  the  respective  power  attenuation  for  each 
pseudolite  was  adjusted  to  maximize  the  SNR  without  getting  signal  cross  correlation 
at  the  receiver.  A  static  position  could  be  computed  once  the  lab  coordinate  system 
was  determined,  the  receiver  time  tag  problem  was  solved,  an  algorithm  to  read  the 
pseudolite  RINEX  data  had  been  constructed,  and  new  data  had  been  taken.  This 
research  worked  toward  a  solution  exclusively  from  pseudoranges,  so  an  accurate 
position  required  an  accurate  pseudorange  derivation  from  the  code  phase  data 
present  in  the  data  file. 

Because  the  time  transmitted  is  not  logged  into  a  navigation  message,  it  has  to 
be  solved  using  the  cycle  counts  inherent  in  the  code  phase  data.  Investigation  of  the 
firmware  code  showed  that  the  resolution  of  the  code  phase  is  2048  samples  per  chip. 
Since  one  code  phase  epoch  {CDP epoch)  has  1023  chips  and  a  rate  of  1 .023e6  chips  per 
second,  there  are  1000  code  phase  epochs  per  second: 
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epoch  Jj,  1 .023e6ChipS  _  ^  ^^^epochs 
IQilZChips  s  s 


(3.1) 


Introducing  the  speed  of  light  (c  =  299792458  m/sec)  results  in  the  code  phase 
length  of  approximately  300km: 


m 

- S - ,  WTiChips  ^  299J92 — - — 

1.023e6?^  CDP,^ 

s 


(3.2) 


A  pseudorange  can  be  computed  as  the  time-received  (TV)  minus  the  time- 
transmitted  (Tf)  in  meters,  times  the  speed  of  light: 

PR  =  iTr-Tt)*c  (3.3) 

The  transmit  time  is  unknown,  but  it  can  be  computed  from  code  phase  since 
the  code  phase  is  sampled  at: 

1 023  *  2048  =  2,095,1 04  =  2,095,1 04,000  (3 .4) 

Chip  CDP,^,,,  s 

Hence,  pseudoranges  can  be  formed  from  the  speed  of  light  times  a  quantity 
containing  the  fractional  part  of  the  received  time  second  minus  the  code  phase 
measurement  (CDP)  scaled  by  the  sampling  resolution  per  second  and  a  factor  that 
represents  the  particular  code  phase  cycle  (A/): 


PR  =  c*(Tr-i- 


CDP 


2,095,104,000 


Samples  1000 


(3.5) 


This  is  the  final  pseudorange  equation  computed  from  code  phase 


measurements.  Because  of  oscillator  run-up,  the  code  phase  cycle  term  (N)  is  a  flag 
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that  changes  the  jigsaw  sealing  effect  of  the  raw  code  phase  data  into  a  continuous 
data  stream.  The  post-proeessed  position  solution  is  then  formed  using  Hofinann- 
Wellenhof  s  linear  model  for  point  positioning  with  pseudoranges  [10].  The  notation 
is  slightly  altered  here  for  consistency  and  simplicity  in  single-receiver  pseudolite 
work,  and  each  term  is  defined  in  Table  2: 

X  —X  Y  —Y  Z  —Z 

R  =p - i - AX — - AT - - AZ  +  c5^-cS  (3.6) 

■'  ^  pj  pj  pj 


J 

Pseudolite  number  (1-4  inNSTL,  1-6  in  simulation) 

Ri 

Measured  code  pseudorange  from  antenna  to  pseudolite 

Pi 

True  geometrie  distance  between  antenna  and  pseudolite 

X,  Yu  Z,- 

NSTL  pseudolite  coordinates 

XY,Z 

Approximate  NSTL  antenna  coordinates 

XX,  AT  AZ 

Corrections  to  approximate  NSTL  antenna  coordinates 

Pseudolite  clock  offset 

5 

Receiver  clock  offset 

c 

Speed  of  Ught  (299792458  m/sec) 

Table  2:  Pseudorange  Point  Positioning  Linear  Model  Term  Deflnitions 

Four  parameters  are  computed  in  the  state  for  the  least  squares  position 
solution:  the  location  (X,  Y,  Z)  of  the  receiver  antenna  with  respect  to  the  lab  origin 
and  the  clock  correction  (cdt).  An  initial  guess  of  (0,  0,  0,  0)  is  used  for  the  antenna 
loeation  and  clock  offset,  and  the  solution  is  calculated  on  an  epoch-by-epoch  basis. 
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Chapter  4:  Code  Phase  Biases 


Initial  investigation  of  the  pseudoranges  calculated  from  the  code  phases 
showed  they  were  of  an  order  of  magnitude  equal  to  that  of  the  speed  of  light,  and 
code  minus  carrier  plots  showed  each  channel  drifting  at  a  different  rate.  This  drift 
could  be  explained  because  the  firmware  code’s  original  frequency  lock  loop  (FLL) 
did  not  directly  track  carrier  phase,  but  instead  tracked  their  derivatives.  A  phase  lock 
loop  (PLL)  was  then  implemented  by  Key  to  directly  track  the  carrier  phases,  and  the 
resulting  code  minus  carrier  plots  showed  similar  drift  rates  for  each  channel  [11]. 

The  first  attempt  at  calculating  a  position  from  code  phase  determined 
pseudoranges  resulted  in  (X,  Y,7)  coordinates  that  quickly  became  undeterminable. 
The  initial  (0,  0,  0, 0)  position  and  clock  estimate  showed  the  problem  became  non¬ 
linear  within  a  few  seconds  because  of  the  close  proximity  (approximately  seven 
meters)  of  the  pseudolites  to  the  receiver  in  comparison  to  a  large  (approximately  300 
kilometers)  pseudorange  measurement.  The  result  was  a  solution  further  from  the 
origin  at  each  successive  epoch.  Trying  a  new  initial  estimate  using  the  speed  of  light 
to  introduce  a  large  order  of  magnitude  in  the  clock  offset  state  parameter  did  not 
have  a  significant  effect  on  linearizing  the  solution. 

In  order  to  keep  the  solution  from  quickly  becoming  undeterminable,  each 
epoch  was  then  referenced  back  to  the  initial  guess  instead  of  the  previous  epoch's 
solution.  These  results  kept  the  solution  within  a  more  "linear"  bound,  but  large  noise 
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factors  still  resulted.  Again,  introducing  the  speed  of  light  into  the  initial  clock 
estimate  did  not  alter  the  solution  at  all. 

The  explanation  for  this  result  was  the  evidence  of  a  bias  inherent  in  each 
pseudolite’s  code  phase.  These  biases  were  also  determined  to  be  somewhat  large. 
Knowing  that  each  pseudolite  was  approximately  seven  meters  from  the  lab  origin,  it 
was  observed  that  the  code  phase  computed  pseudorange  position  solution  was 
varying  in  the  fourth  order  of  magnitude  because  the  oscillator  to  which  all  of  the 
pseudolites  were  referenced  did  not  use  the  fundamental  GPS  frequency  (1 .023  MHz) 
inherent  in  the  receiver.  Instead,  the  reference  oscillator  used  20.46  MHz.  This 
meant  that  each  pseudolite  had  the  opportunity  to  begin  its  code  phase  data  stream  at 
one  of  twenty  different  times  within  the  transmit  time  second.  Hence,  these  biases 
had  to  be  determined  and  removed  before  these  pseudoranges  could  be  of  any  use. 

The  first  attempt  at  removing  the  biases  was  unsuccessful.  A  second  order 
polynomial  was  closely  fit  to  a  common  pseudolite’s  code  phase  and  was  removed 
from  all  pseudolites.  It  was  originally  thought  that  the  remaining  information  was  the 
resulting  code  phase  for  each  respective  pseudolite.  However,  once  these  code  phases 
were  transformed  into  pseudoranges  (discussed  in  Chapter  3)  and  later  smoothed 
using  carrier  phases  (later  discussed  in  Chapter  6),  the  expected  white  Gaussian  noise 
was  not  the  only  error  that  remained.  Consequently,  removing  a  common  polynomial 
from  each  code  phase  did  not  accoimt  for  the  fact  that  each  bias  was  unique  to  its 
corresponding  pseudolite.  As  a  result,  a  Kalman  filter  approach  was  implemented  to 

remove  each  bias  before  the  receiver’s  position  was  computed. 
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Chapter  5:  Bias  Kalman  Filter 


The  overall  goal  of  this  work  was  to  achieve  a  single  post  processed  position 
solution  within  five  meters  of  the  true  position  due  to  the  Gaussian  white  noise  on  the 
computed  pseudoranges.  However,  a  Kalman  filter  approach  to  resolving  the  code 
phase  biases  was  chosen  as  opposed  to  a  batch  algorithm  so  that  future  work  can 
eventually  solve  for  a  real-time  solution  and  eliminate  the  post-processing 
requirement.  Based  on  the  results  presented  in  Chapter  4,  these  biases  were  assumed 
to  resemble  a  time  invariant  offset;  each  one  unique  and  inherent  to  a  pseudolite’s 
code  phase.  At  this  point  in  the  research,  four  pseudolites  remained  fixed  in  the 
NSTL  because  two  of  the  original  six  were  removed  for  closer  study  at  CSR. 
Therefore,  it  was  possible  to  estimate  three  of  the  four  biases  and  dedicate  one  of  the 
code  phases  as  the  reference  datum.  In  other  words,  one  bias  was  assumed  to  be  zero 
and  the  other  three  were  estimated  using  it  as  a  reference.  By  including  these  biases, 
the  pseudorange  measurement  equation  (to  each  pseudolite  J)  took  the  form: 

PRj  =  pj  +cdt  +aj+T]  (5.1) 

In  this  equation,  the  clock  offsets  fi’om  the  receiver  and  reference  oscillator 
were  combined  into  one  bias  (cdt)  common  to  all  four  pseudoranges.  The  first  code 
phase  bias  (aj)  equaled  zero,  and  there  was  a  unique  geometric  true  range  (pj) 
between  the  receiver  and  each  pseudolite.  The  last  parameter  (rj)  represented  the 
remaining  unknown  noise  on  the  measurement  that  would  theoretically  carry  through 
to  the  position  solution. 
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Before  building  a  Kalman  filter  that  solved  for  unknown  biases  in  NSTL  data, 
it  was  beneficial  to  filter  simulated  data  first  and  ensure  that  known  biases  could  be 
resolved.  Using  a  pseudolite  data  simulator  (PSD_SIM.EXE)  written  by  Key, 
pseudorange,  carrier  phase,  smoothed  pseudorange,  and  Doppler  data  streams  subject 
to  the  same  measurable  characteristics  observed  in  the  NSTL  were  created  [12].  The 
simulator  returned  RINEX  data  files  similar  to  NSTL  tests.  The  simulated 
pseudorange  (PR)  information  was  primarily  usefixl  to  test  the  filtering,  but  the 
simulated  carrier  phase  (<D)  information  was  also  used  to  test  the  smoothing  algorithm 
that  will  be  discussed  in  Chapter  6.  The  equations  used  to  produce  the  simulated  data 
were  very  similar  to  those  modeled  in  the  NSTL: 


PRj  =  pj  +  c{dT  -  dt)  +  +  // 

=  Pj  +  c{dT  -  dt)  +  XN j  +  V 


(5.2) 


This  simulator  allowed  the  user  to  define  the  output  pseudorange  and  carrier 
phase  noise  (rj  and  o  respectively).  It  also  allowed  the  user  to  toggle  whether  the 
receiver  position  was  offset  fi'om  the  lab  origin  or  in  motion,  and  whether  the 
measurements  included  line  biases  (between  the  reference  oscillator  and  the 
pseudolites),  code  and  carrier  phase  biases  (a  and  respectively),  or  reference  and 
receiver  clock  offset  and  drift  {dT  and  dt  respectively).  The  simulated  clock  offset 
values  were  formed  fi-om  a  randomly  generated  bias  and  fi-equency,  and  the  simulated 
true  geometric  ranges  (p)  were  consistent  with  those  measured  in  the  NSTL 
theodolite  survey.  Constant  known  code  phase  biases  were  initially  added  to  “clean” 
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simulated  pseudoranges  (with  minimal  biases  and  noise),  and  the  Kalman  filter 
resolved  those  biases  perfectly.  The  simulated  biases  and  noises  were  then  added  one 
at  a  time  to  determine  how  the  a  priori  covariance  and  normal  noise  distribution 
matrices  should  be  adjusted  as  defined  in  the  Kalman  filter  algorithm  below. 

Four  pseudorange  equations  were  modeled  using  NSTL  data  while  six 
equations  were  modeled  using  simulated  data.  There  was  no  significant  difference  in 
the  filtering  by  losing  two  pseudohtes.  However,  estimating  the  three  position 
parameters  (X,  Y,  Z),  the  clock  offset,  and  the  three  code  phase  biases  at  one  epoch 
resulted  in  seven  state  parameters.  A  singularity  problem  occurred  because  there 
were  four  equations  and  seven  unknowns.  Adding  more  pseudohtes  did  not  help  this 
problem  because  each  additional  pseudoUte  would  add  another  bias  along  with  its 
measurement.  The  singularity  occurred  while  computing  the  (HH^)''  portion  of  the 
gain  matrbc.  Because  a  time  varying  parameter  such  as  motion  (velocity  or 
acceleration)  was  not  modeled  in  this  scenario,  the  observation  matrbc  did  not  change 
fi'om  epoch  to  epoch.  The  HH^  computation  was  not  fuU  rank,  and  therefore  its 
inverse  did  not  exist  when  there  were  more  unknowns  than  equations.  Also,  each 
measurement  was  weighted  the  same,  and  the  state  transition  matrbc  used  in  the  filter 
did  not  update  any  type  of  dynamics  over  time.  Both  of  these  parameters  were 
therefore  modeled  as  the  identity  matrix  throughout  the  filtering. 

Because  of  the  singularity  issue,  the  final  solution  to  this  problem  required 
two  Kalman  filters.  The  first  was  considered  a  cahbration  filter  used  to  determine  the 

constant  code  phase  biases  and  clock  ofiset  during  a  period  over  a  known  position. 
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The  second  filter  was  used  after  the  calibration  period  to  solve  for  the  receiver’s 
position  once  the  code  phase  biases  were  known.  In  contrast  to  the  receiver  clock 
offset  and  drift  evident  in  the  solution,  each  code  phase  bias  was  different  for  each 
pseudolite.  Hence  for  each  epoch,  the  clock  bias  term  was  subtracted  from  each 
measurement  in  addition  to  the  unique  code  phase  bias  for  that  pseudolite. 

Both  filters  were  derived  using  the  algorithms  written  by  Gelb  and  Tapley 
[13]  [14].  Constructed  using  Matlab  scripts,  the  source  code  for  the  NSTL  data 
solution  can  be  found  in  the  Appendix.  The  simulated  data  position  algorithm  was 
only  slightly  different  than  the  NSTL  data  algorithm  and  was  therefore  not  listed  here. 
It  only  varied  in  that  pseudoranges  were  not  computed  from  code  phases  since  they 
were  directly  provided  from  the  simulated  RINEX  output  file,  and  the  number  of 
measurements  was  increased  from  four  to  six.  The  entire  NSTL  data  position 
algorithm  is  described  as  follows. 

1 .  Read  the  raw  code  phases  from  RINEX  file. 

2.  Compute  the  raw  (no  repeating)  pseudoranges  (PR’s)  from  the  code  phases. 

3.  Initialize  the  calibration  nominal  sate.  Seven  meters  is  the  approximate  range 
from  the  receiver  to  each  pseudolite,  cdt  represents  the  entire  clock  bias 
(offset  and  drift),  and  a  represents  the  respective  code  phase  biases: 


cdt  PR^  -  7 
a2  PR2  -1  -cdt 
PR2  -  7  -  cdt 
-1  -cdt 


(5.3) 
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4.  Initialize  the  calibration  a  priori  covariance  matrix  (Pq),  normal  noise 


distribution  matrix  (5),  state  transition  matrix  (O),  and  weight  matrix  (W): 


'10 

0 

0 

o' 

'1 

0 

0 

0  ■ 

0 

2 

0 

0 

s= 

0 

0.05 

0 

0 

0 

0 

2 

0 

0 

0 

0.05 

0 

0 

0 

0 

2_ 

0 

0 

0 

0.05_ 

'l 

0 

0 

o' 

"1 

0 

0 

o' 

0)  = 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

w  = 

1 

0 

0 

0 

1 

0 

0 

0 

0 

1_ 

0 

0 

0 

1_ 

5.  Loop  over  all  calibration  epochs  of  data. 

a.  Read  epoch. 

b.  Propagate  the  covariance  matrix: 

P*  +S 


(5.4) 


(5.5) 


c.  Calculate  the  observation-state  (H)  and  observation  (y)  matrices  where 


p  is  the  respective  (X,  Y,  Z)  norm  from  the  receiver  to  each  pseudolite: 


'l 

0 

0 

o' 

PR,  -/?,  -cdt 

1 

1 

0 

0 

PR2  -  P2-  cdt 

1 

0 

1 

0 

y  = 

PR^  -  P3  -  cdt 

1 

0 

0 

1 

PR^  -  P4  -  cdt 

d.  Calculate  the  gain  matrix: 

K  =  P,H^[hP,H^  (5.7) 

e.  Calculate  the  new  estimate: 


A 


x  =  Ky 
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(5.8) 


f.  Update  the  state  and  covarianee  matrix: 


Z  =  J^  +  jc  P^_,=P^-KHP,^ 


(5.9) 


g.  Read  the  next  epoeh. 

6.  End  the  calibration  loop. 

7.  Initialize  the  position  nominal  sate: 


'X' 

■  0.0  * 

Y 

z= 

0.0 

Z 

-0.2_ 

(5.10) 


8.  Initialize  the  position  a  priori  eovariance  matrix  (Pq),  normal  noise 

distribution  matrk  (5),  state  transition  matrk  (O),  and  weight  matrix  (W): 


2 

0 

0' 

'0.05 

0 

0 

Po  = 

0 

2 

0 

s  = 

0 

0.05  0 

0 

0 

2_ 

0 

0 

0.05  J 

'1 

0 

O' 

'1 

0 

o' 

<D  = 

0 

1 

0 

w  = 

= 

0 

1 

0 

0 

0 

1 

0 

0 

1 

9.  Loop  over  all  position  epoehs  of  data. 

a.  Read  epoch. 

b.  Propagate  the  covariance  matrix: 

P*  +5  (5.12) 

c.  Calculate  the  observation-state  (H)  and  observation  (y)  matrices: 
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X'  -X  r‘ -y  Z'  -z' 

Pj  P\  P\ 

-X  -y  -z  r  PR,-p^-cdt 

P2  Pi  Pi  _  PR-I  ~  Pi  ~  ~  ®2 

X^  -X  -y  Z" -z  PR^-p^-cdt-a^ 

Pi  Pi  Pi  PR4  -  P4  -  cdt  - 

X'*  -X  Y^ -y  Z^ -z 

Pi,  Pi  Pi 

d.  Calculate  the  gain  matrix: 

K  =  P^H^[hP^H^  +W-^Y  (5.14) 

e.  Calculate  the  new  estimate: 

x  =  Ky  (5.15) 

f.  Update  the  state  and  covariance  matrix: 

X  =  X  +  i  P,_,=P,-KHP,  (5.16) 

g.  Read  the  next  epoch. 

10.  End  the  position  loop. 

Using  this  algorithm  to  estimate  a  position  solution  proved  to  be  very  accurate 
in  simulation.  Using  thirty  minutes  of  simulated  data  with  two  meters  of  noise  on  the 
pseudoranges,  a  random  clock  offset,  random  code  phase  biases,  and  a  receiver 
location  (0.029  m,  -0.016  m,  -0.172  m)  slightly  offset  form  the  origin,  the  following 
results  were  obtained.  Figure  4  illustrates  the  calibration  clock  bias  and  code  phase 
biases  for  the  sk  simulated  pseudolites. 
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Figure  4:  Simulated  Clock  Bias  and  Code  Phase  Biases 


The  magnitudes  of  these  biases  resembled  what  were  expected  from  NSTL 
data.  In  each  case,  the  filter  did  not  converge  to  better  than  two  meters  of  noise  on 
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each  bias.  It  was  therefore  assumed  that  the  bias  value  obtained  at  the  end  of  the  filter 
was  accurate  to  within  the  noise  inherent  in  the  pseudoranges.  Figure  5  illustrates  the 
simulated  position  solution  estimate  after  subtracting  the  biases  in  Figure  4. 

X  Estimate 


Figure  5:  Simulated  Position  Solution 


These  results  showed  approximately  two  meters  of  noise  and  the 

characteristics  of  the  original  receiver  location  (0.029  m,  -0.016  m,  -0.172  m). 

Consequently,  the  calibration-position  algorithm  was  validated  for  simulated  data  that 

closely  resembled  what  was  expected  in  the  NSTL. 

Figures  6  and  7  illustrate  the  unknown  biases  and  position  solution  estimated 

fi-om  NSTL  data.  This  data  file  (NSTL1930.000)  contained  an  hour  of  data  taken  by 

Lightsey  [15].  The  receiver’s  location  was  approximately  (0.0  m,  0.0  m,  -0.2  m). 

The  negative  Z-coordinate  described  a  slight  elevation  above  the  origin  because  the 

lab  coordinate  system  was  defined  using  positive  values  in  the  “down”  direction. 
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Figure  6:  NSTL  Clock  Bias  and  Code  Phase  Biases 
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Figure  7:  NSTL  Position  Solution 
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Unfortunately,  it  was  obvious  that  the  NSTL  position  solution  was  not  nearly 
as  aecurate  as  that  formed  in  simulation.  The  solution  wandered  on  the  order  of  30 
meters  in  each  direction.  Possible  reasons  for  such  errors  could  have  been  from  the 
assumption  that  the  code  phase  biases  were  constant.  Studying  the  code  phase  bias 
estimates  in  the  NSTL  data  (Figure  6)  showed  a  maximum  root  sum  squared  (RSS) 
variation  in  of  approximately  75  meters  around  an  absolute  mean  throughout  the 

data  file.  However,  the  standard  deviations  (the  square  roots  of  the  variances  in  the 
covariance  matrbc)  on  all  four  bias  estimates  {cdt,  ai,  as,  and  a4)  were  approximately 
55  centimeters  over  aU  of  the  epochs.  This  meant  that  Figure  6  possibly  illustrates  the 
true  bias  at  each  epoch  and  not  just  random  noise.  Thus,  the  biases  may  have  been 
changing  over  time,  or  multipath  could  have  been  affecting  the  resulting  bias 
estimates.  The  error  seen  in  the  position  solution  could  have  therefore  directly 
transferred  from  the  errors  in  the  bias  estimates.  Because  this  solution  was  not  very 
useful  in  its  current  form,  a  smoothing  technique  using  carrier  phases  was 
implemented  after  the  filtering  in  an  attempt  to  achieve  a  more  precise  position. 
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Chapter  6:  Carrier  Phase  Smoothed  Code 


Lachapelle’s  pseudorange  smoothing  using  carrier  phases  (outlined  below 
from  Hofi&nan-WellenhoQ  relies  on  the  increased  accuracy  provided  by  the  carrier 
phase  data  over  the  noisy  pseudoranges  [10].  In  order  to  handle  the  large  errors  in 
pseudorange  data,  the  pseudoranges  (PR)  are  weighted  more  heavily  during  early 
epochs  and  the  carrier  phases  (d))  are  weighted  more  heavily  during  later  epochs.  At 
the  first  epoch,  the  weight  (w)  equals  one,  and  it  is  decreased  for  each  successive 
epoch  by  an  appropriate  fector  determined  by  the  length  of  the  data  file  until  the 
dependence  on  the  pseudoranges  is  negligible  or  zero.  Although  cycle  shps  have  not 
been  found  to  readily  occur  in  NSTL  data  files,  should  one  occur,  the  weight  would 
be  reset  to  one  and  slowly  include  the  carrier  phases  again.  This  smoothed 
pseudorange  equation  can  be  computed  as; 

PR^sm  =  wPR^  +  (l  -  w\PR^_^sm  +  )  (6.1) 

For  the  work  done  here,  full  weight  was  given  to  the  carrier  phases  after  half 
of  the  data  was  processed  in  order  to  obtain  precision  in  the  position  information 
provided  by  the  pseudoranges.  This  procedure  was  coded  using  a  Matlab  script  in 
conjunction  with  a  least  squares  position  solution  after  the  pseudoranges  were 
computed  from  the  code  phases  (described  in  Chapter  5)  and  the  code  phase  biases 
were  removed  (described  in  Chapter  6).  The  NSTL  data  pseudorange  smoothing 
code  can  be  seen  in  the  Appendix.  It  should  be  noted  that  in  kinematic  or  real-time 
applications,  more  weight  could  be  given  to  the  carrier  phases  at  an  earlier  epoch. 
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Figure  8  shows  how  this  smoothing  technique  worked  very  well  in  simulation. 
Using  the  same  simulated  data  file  that  produced  the  biases  and  position  estimate  in 
Figures  4  and  5,  this  technique  resulted  in  centimeter-level  error  fi-om  the  actual 
simulated  receiver  location  (0.029  m,  -0.016  m,  -0.172  m). 

Sim  XPos,  Smoothed  A\erage:  0.41875m.  Forward  Smoothed,  Weight;  1->0,  Step:  0.001 


0  200  400  600  800  1000  1200  1400  1600  1800 


Sim  Y  Pos,  Smoothed  Average:  -0.41786m 


0  200  400  600  800  1000  1200  1400  1600  1800 

Sim  Z  Pos,  Smoothed  Average:  0.14359m 


0  200  400  600  800  1000  1200  1400  1600  1800 

Epoch  Second  (#  of  observations) 

Figure  8:  Simulated  Smoothed  Position  Solution 

In  NSTL  data,  this  smoothing  technique  also  proved  to  be  somewhat  useful, 
but  not  ideal.  Figure  9  shows  smoothed  pseudoranges  fi-om  data  file  NSTL1930.000 
after  the  biases  in  Figure  6  were  removed.  The  resulting  error  was  approximately  five 
meters  fi-om  the  (0.0  m,  0.0  m,  -0.2  m)  actual  receiver  location. 
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Lab  XPos,  Smoothed  Average:  -0.14489m,  Forward  Smoothed,  Weight:  1->0,  Step:  0.0005 


Figure  9:  NSTL  Smoothed  Position  Solution 


While  the  NSTL  smoothed  solution  was  better  than  the  position  estimate 
provided  by  the  Kalman  fQter,  it  was  still  not  to  the  expected  centimeter-level 
accuracy  that  was  achieved  in  simulation.  The  combination  of  filtering  and 
smoothing  the  pseudoranges  was  validated  in  simulation,  but  there  could  still  have 
been  large  multipath  issues  or  time-varying  code  phase  biases  that  caused  the  NSTL 
position  estimate  to  be  less  accurate. 
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Chapter?:  Conclusion 


Achieving  a  static  position  solution  from  an  indoor  pseudolite  system  such  as 
the  one  studied  here  in  the  NSTL  is  more  complicated  than  normal  GPS  in  certain 
aspects.  Because  the  code  phases  were  not  perfeetly  synchronized,  an  extra  bias 
existed  in  the  pseudorange  from  each  pseudolite.  A  two-step  Kalman  filter  and 
carrier  phase  smoothing  was  required  to  resolve  these  biases.  The  first  filter  was  used 
during  a  calibration  period  over  a  known  location  so  the  constant  biases  could  be 
determined.  The  second  filter  solved  for  a  position  once  the  biases  were  known. 

In  simulation,  this  method  worked  extremely  well.  The  filter  constructed  here 
was  successful  in  determining  the  biases  within  a  few  eentimeters  in  simulated  data. 
The  simulated  data  closely  resembled  what  was  expeeted  in  real  NSTL  data  beeause 
in  addition  to  the  code  phase  biases,  it  had  pseudorange  noise  and  clock  biases  similar 
in  magnitude  to  those  seen  in  the  NSTL.  Although  the  simulation  was  designed  to 
model  a  certain  interpretation  of  the  measurement  environment,  it  does  not  appear  to 
be  perfectly  consistent  with  what  is  observed  in  the  NSTL. 

The  jump  from  simulated  to  NSTL  data  approximately  resulted  in  a  30-meter 
increase  in  bias  and  position  error.  Smoothing  using  the  accuracy  provided  by  the 
carrier  phases  was  then  implemented.  When  both  the  simulated  and  NSTL  position 
estimates  were  smoothed  in  this  manner,  the  simulated  data  found  the  receiver 
location  within  a  few  centimeters  while  the  NSTL  solution  was  accurate  to  five 
meters.  Table  3  eompares  the  root  sum  squared  (RSS)  range  of  errors  in  the  filtered 
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position  solutions  and  the  root  mean  squared  (RMS)  errors  in  the  final  smoothed 
position  solutions  fi-om  the  receiver’s  true  location  for  both  simulated  and  NSTL  data. 


Simulation 

NSTL 

RSS  on 

X 

2.3 

32.0 

Calibration/Position 

Y 

2.2 

31.0 

Estimate 

Z 

2.1 

32.0 

RMS  on 

X 

+  0.39 

-0.14 

Final  Smoothed 

Y 

-0.40 

+  2.09 

Solution 

Z 

+  0.32 

+  4.50 

Table  3:  Position  Error  Summary  (Meters) 

While  the  NSTL  position  solution  fi-om  this  method  was  not  perfect,  it  was  not 
altogether  insignificant.  The  noise  and  bias  characteristics  in  the  lab  data  were  not 
perfectly  understood,  and  therefore  may  not  have  been  perfectly  estimated.  More 
specifically,  it  could  have  been  possible  that  the  code  phase  biases  were  not  constant. 
They  could  have  varied,  as  does  multipath  with  time.  Multipath  itself  coxild  have 
introduced  a  large  portion  of  the  30-meter  error  that  forced  the  position  to  change 
with  time  because  the  NSTL  was  not  a  controlled  environment.  Large  metal  doors 
and  equipment,  the  operation  of  nearby  electronic  equipment,  and  people’s 
movements  within  the  lab  could  have  all  magnified  the  effects  of  multipath  while  the 
tests  were  conducted. 

Moreover,  the  hardware  or  the  post-processing  algorithm  may  have  been 
inaccurate.  If  the  pseudolites  and  the  oscillator  were  not  exactly  synchronized  in  the 
manner  described  in  Chapters  2  and  3,  time  varying  biases  may  have  been  produced 
by  the  hardware.  Otherwise,  the  understanding  of  the  firmware  code  that  resulted  in 
the  code  phase  to  pseudorange  derivation  in  Chapter  3  may  have  been  incorrect. 
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Depending  on  the  magnitude  of  the  problem,  erroneous  pseudoranges  could  have 
resulted  in  additional,  possibly  time-varying,  biases  that  were  not  properly  modeled  in 
the  filtering. 

In  any  event,  it  is  recommended  that  the  next  step  is  to  validate  the  variability 
in  the  code  phase  biases.  As  an  example,  motion  could  be  apphed  to  the  receiver 
after  the  calibration  period  to  determine  if  that  motion  can  be  modeled  using  the 
constant  biases  resolved  by  this  filter.  If  the  five-meter  noise  is  still  apparent  on  the 
smoothed  solution  but  the  motion  can  be  deciphered,  then  the  noisy  results  seen  here 
may  very  well  be  the  best  available  solution  under  the  NSTL  environment. 

If  a  motion  test  does  not  isolate  the  situation  as  a  time- varying  bias  problem, 
then  the  answer  may  be  that  the  code  phases  must  be  properly  synchronized  in  order 
to  get  a  static  position  solution  using  only  the  pseudoranges.  In  reality,  this  would  be 
a  significant  recommendation  to  negate  the  errors  inherent  in  the  filtering. 

Regardless,  the  work  accomplished  here  did  achieve  a  five-meter  solution  in 
an  extremely  noisy  environment.  In  space,  an  entirely  different  environment  will  be 
encoimtered,  and  new  problems  will  most  likely  appear.  Nonetheless,  this  work  has 
provided  valuable  insight  into  pseudolites  and  is  a  stepping-stone  to  eventual 
kinematic  pseudolite  applications  in  space. 
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Appendix 


%Scott  Weyermiiller 

%Code  Phase  Bias  Kalman  Filter  During  Calibration  Period 
%Open  observation  file  using  startpseud.m  and  readcdp.m  to  read  PR's 
%prior  to  running  this  script 
%Last  Updated  3  Aug  2000 

%  ‘  . . . . . 

clear  X_0  X  xhat  Pbar  Pkminl  Pk  P  S  H  y  Phi  W  K 
clear  outH  outy  outK  outxhat  outX  outP  Obs  beg  obstime 


%Set  constants 
c=299792458.0; 
pseud-size(sv,2); 
obstime=size(time,  1 ); 


%speed  of  light  in  m/s 

%pseud  is  #  of  observed  pseudolites  (4) 

%obstime  is  #  of  observation  times  (3379) 


%Read  NSTL  pseudolite  eoordinates 
plabl93; 


%Initialize  all  matrix  dimensions 
X_0=zeros(4,l); 

X=zeros(4,l); 
xhat=zeros(4,l); 
Pbar=zeros(4,4); 

Pkminl  =zeros(4, 4); 
Pk=zeros(4,4); 

P=zeros(4,4); 

S=zeros(4,4); 

H=zeros(4,4); 
y=zeros(4,l); 

K=zeros(4,4); 

%Set  initial  value 
dti=PR(l,l)-7; 
a2i=PR(l,2)-7-dti; 
a3i=PR(l,3)-7-dti; 
a4i=PR(l,4)-7-dti; 

X_0=[dti,  a2i,  a3i,  a4i]; 


%lnitial  guess  of  state  values 
%New  state  after  each  epoch 
%State  update  after  each  epoch 
%Apriori  covariance  matrix 
%Covariance  matrix  of  last  epoch 
%Covariance  matrix  of  current  epoch 
%Covariance  matrix  holder 
%Normal  noise  distribution  matrix 
%Model  matrix 

%Observed-Computed  equation 
%Gain  matrix 


of  nominal  X  (dti,  a2i,  a3i,  a4i) 


%Fill  apriori  covariance  matrix 
Pbar(l,l)=10; 
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Pbar(2,2)=2; 

Pbar(3,3)=2; 

Pbar(4,4)=2; 

Pkminl=Pbar;  %Last  epoch  covariance  equals  apriori  for  first  computation 

%Fill  S  matrix  to  open  noise  search  parameters  in  covariance  matrix 

s(i,i)=i; 

S(2,2)=.05; 

S(3,3)=.05; 

S(4,4)=.05; 

%Fill  estimate,  state  transition  and  weight  matrices 
X(1 :4)=X_0;  %First  epoch  uses  initial  guess  of  state 

Phi=eye(4);  %State  transition  matrix  equals  identity 

W=eye(4);  %Weight  matrix  equals  identity 

calpos=[0.0  0.0  -0.2]; 

%Loop  over  all  observation  times 
for  oloop=l:obstime 

%Propogate  covariance  matrix  to  current  epoch 
Pk=Phi*Pkminl  *Phi'+S; 

%Begin  available  pseudolite  loop 
for  ploop=l  rpseud 

%Calculate  range  from  last  receiver  position  (estimate) 
posdif!=pl(ploop,  1  ;3)-calpos; 
range=norm(posdiff); 

%Fill  H  matrix  with  clock  offset  and  bias  values 

H(ploop,l)=l ;  %Use  1  here  instead  of  c  to  avoid  singularity  issues  later 

H(2,2)=l; 

H(3,3)=l; 

H(4,4)=l; 

%Fill  y  vector  (observed-computed)  for  all  pseudorange  values 
if  ploop==l  %First  PR  bias  equals  zero,  so  don't  include  here 

y(ploop)=PR(oloop,ploop)-range-X(  1 ) ; 
else 

y(ploop)=PR(oloop,ploop)-range-X(  1  )-X(ploop); 
end 
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%End  available  pseudolite  loop 
end 

%Calculate  gain  matrix 
K=Pk*H'*inv(H*Pk*H'+inv(W)); 

%Calculate  new  estimate  and  update  estimated  state 
xhat=K*y; 

X=X+xhat; 

%Update  covariance  matrix 
P=Pk-K*H*Pk; 

Pkminl=P; 

%Store  output  values 

outH( : , :  ,oloop)=H; 

outy(:,:,oloop)=y; 

outK( : , :  ,oloop)=K; 

outxhat(:,:,oloop)=xhat; 

outX(:,:,oloop)=X; 

outP(:,:,oloop)=P; 

Obs(oloop)=oloop; 

%End  observation  time  loop 
end 

%Manipulate  stored  output  values  for  plots 

update=squeeze(outxhat);  %Holds  bias  updates  for  each  epoch 

uClkrec=update(l 

ua2=update(2,:); 

ua3=update(3,:); 

ua4=update(4,;); 

posn=squeeze(outX);  %Holds  updated  state  for  each  epoch 

Clkrec=posn(l,:); 

a2=posn(2,:); 

a3=posn(3,:); 

a4=posn(4,:); 

Xfinal=X; 

%Plot  results 
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%Scott  Weyermuller 

%X,Y,Z  Position  Kalman  Filter  After  Calibration  Period 

%Open  observation  file  using  startpseud.m,  form  PR’s  using  readcdp.m,  and  find 
%code  phase  biases  using  kalmanbias.m  prior  to  running  this  script 
%Last  Updated  3  Aug  2000 


clkrem=outX; 

clear  X_0  X  xhat  Pbar  Pkminl  Pk  P  S  H  y  Phi  W  K  outH  outy  outK  outxhat  outX 
outP  Obs  beg  obstime  posdiff  range 


Clrem(:,l)=PR(;,l); 

Clrem(:,2)=PR(:,2)-Xfinal(2); 

Clrem(;,3)=PR(:,3)-Xfinal(3); 

Clrem(:,4)=PR(:,4)-Xfinal(4); 

%Set  constants 
c=299792458.0; 
pseud=size(sv,2); 
obstime=size(time,  1 ); 


%speed  of  light  in  m/s 

%pseud  is  #  of  observed  pseudolites  (4) 

%obstime  is  #  of  observation  times  (3379) 


%Read  NSTL  pseudolite  coordinates 
plabl93; 

%Initialize  all  matrix  dimensions 
X_0=zeros(3,l); 

X=zeros(3,l); 

xhat=zeros(3,l); 

Pbar=zeros(3,3); 

Pkminl=zeros(3,3); 

Pk=zeros(3,3); 

P=zeros(3,3); 

S=zeros(3,3); 

H=zeros(4,3); 

y=zeros(4,l); 

K=zeros(3,4); 


%Tnitial  guess  of  state  values 
%New  state  after  each  epoch 
%State  update  after  each  epoch 
%Apriori  covariance  matrix 
%Covariance  matrix  of  last  epoch 
%Covariance  matrix  of  current  epoch 
%Co variance  matrix  holder 
%Normal  noise  distribution  matrix 
%Model  matrix 

%Observed-Computed  equation 
%Gain  matrix 


%Set  initial  value  of  nominal  X  (xi,  yi,  zi) 

xi=0.0; 

yi=0.0; 

zi=-0.2; 

X_0=[xi,  yi,  zi]; 
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%Fill  apriori  covariance  matrix 
Pbar(l,l)=2; 

Pbar(2,2)=2; 

Pbar(3,3)=2; 

Pkminl=Pbar;  %Last  epoch  covariance  equals  apriori  for  first  computation 

%Fill  S  matrix  to  open  noise  search  parameters  in  covariance  matrix 
S(l,l)=05; 

S(2,2)=.05; 

S(3,3)=.05; 

%Fill  estimate,  state  transition  and  weight  matrices 
X(1 :3)=X_0;  %First  epoch  uses  initial  guess  of  state 

Phi=eye(3);  %State  transition  matrix  equals  identity 

W=eye(4);  %Weight  matrix  equals  identity 

%Loop  over  all  observation  times 
for  oloop=l:obstime 

%Propogate  covariance  matrix  to  current  epoch 
Pk=Phi*Pkminl  *Phi'+S; 

%Begin  available  pseudolite  loop 
for  ploop=l  :pseud 

%Calculate  range  from  last  receiver  position  (estimate) 
posdiff=pl(ploop,l  :3)-X(l  :3)'; 
range=norm(posdiff); 

%Fill  FI  matrix  with  range  values 
H(ploop,  1  )=-(pl(ploop,  1  )-X(  1  ))/range; 

H(ploop,2)=-(pl(ploop,2)-X(2))/range; 

H(ploop,3)=-(pl(ploop,3)-X(3))/range; 

%Fill  y  vector  (observed-computed)  for  all  pseudorange  values 
y(ploop)=C  1  rem(oloop,ploop)-range-clkrem(l ,  1  ,oloop); 

%End  available  pseudolite  loop 
end 

%Calculate  gain  matrix 
K=Pk*H'*inv(H*Pk*H'+inv(W)); 
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%Calculate  new  estimate  and  update  estimated  state 
xhat=K*y; 

X=X+xhat; 

%Update  covariance  matrix 
P=Pk-K*H*Pk; 

Pkminl=P; 

%Store  output  values 

outH(:,:,oloop)=H; 

outy(:,:,oloop)=y; 

outK(:,:,oloop)=K; 

outxhat(:,:,oloop)=xhat; 

outX(:,:,oloop)=X; 

outP(:,:,oloop)=P; 

Obs(oloop)=oloop; 

%End  observation  time  loop 
end 

%Manipulate  stored  output  values  for  plots 

update=squeeze(outxhat);  %Holds  updates  for  each  epoch 

uposx=update(  1 , 

uposy=update(2,:); 

uposz=update(3,:); 

posn=squeeze(outX);  %Holds  updated  state  for  each  epoch 

posx=posn(l,:); 

posy=posn(2,:); 

posz=posn(3,:); 

Xfmalpos=X; 

%Plot  results 
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%Scott  Weyermuller 

%GPS  Pseudolite  Point  Positioning  Algorithm  Using  Lab  Carrier-Smoothed  PR’s 
%Run  startpseud.m  and  readcdp.m  prior  to  running  this  script  to  read  PR's 
%Specify  which  plab.m  PL  coordinates  to  use  and  which  type  of  position 
%(C1,  PR,  smPR,  Clsm,  etc.)  to  calculate 
%Last  Updated  3  Aug  2000 

clear  A 1 Q  dx; 

file='nstll  931.000'; 
c=299792458.0; 
n=size(sv,2); 
m=size(time,l); 
w  =  1; 

wvar  =  0.0005; 
wmin  =  0.0; 

global  file  m  n 

%Call  smoothing  function  and  read  PL  lab  coordinates 
[Clsm,  avg]  =  smooth(Clrem,  LI,  w,  wvar,  wmin); 
plabl93; 

%lnitial  receiver  position  and  clock  guess 
pos0=[xi,  yi,  zi,  dti]; 

%Begin  observation  time  loop 
for  iloop=l  :m 

%Begin  available  sat  loop 
for  j=l;n 

%Calculate  range  from  last  receiver  position 
pos0diff=pl(j,l  :3)-pos0(l  :3); 
range=norm(pos0difi); 

%Fill  A  matrix  for  all  sat's  (let  c==l) 

A(j,  1  )=-(pl(j.  1  )-pos0(  1  ))/range; 

A(j,2)=-(plO,2)-posO(2))/range; 

A0,3)=-(pia,3)-pos0(3))/range; 

A(j,4)=l; 


%speed  of  light  in  m/s 

%n  is  #  of  observed  pseudolites  (4) 

%m  is  #  of  observation  times  (3379) 
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%Fill  1  vector  for  all  pi's  when  bias  is  estimated 
l(j)=C  1  sm(iloop,j)-range-clkrem(  1 , 1  ,iloop); 

%End  available  sat  loop 
end 

%Fill  cofactor  matrix 
Q=inv(A'*A); 

%Compute  gdop  and  dx 

gdop=sqrt(Q(  1 , 1  )+Q(2,2)+Q(3,3)+Q(4,4)); 

dx=Q*A'*r; 

%Update  receiver  position  for  successive  iterations 

pos(  1  )=posO(  1  )+dx(  1 ); 

pos(2)=^os0(2)+dx(2); 

pos(3)=pos0(3)+dx(3); 

pos(4)=pos0(4)+dx(4); 

%Store  values  for  plots 
Posx(iloop)=pos(  1 ) ; 

Posy(iloop)=pos(2); 

Posz(iloop)=pos(3 ) ; 

Clkrec(iloop)=pos(4); 

GDOP(iloop)=gdop; 

DOPx(iloop)=Q(  1,1); 

DOPy(Uoop)=Q(2,2); 

DOPz(iloop)=Q(3,3); 

DOPt(iloop)=Q(4,4); 

Obs(iloop)=iloop; 

%End  observation  time  loop 
end 

nix=num2str(mean(Posx(avg:m))); 

my=nuin2str(mean(Posy(avg:m))); 

mz=nuni2str(mean(Posz(avg:m))); 

w=num2str(w); 

wmin=num2str(wniin) ; 

wvar=num2str(wvar); 

%Plot  values  vs  #  of  observations 
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%Scott  Weyermuller 

%PR  Smoothing  Function  from  Carrier  Phases 

%Run  startpseud.m  and  readcdp.m  prior  to  running  this  script  to  read  PR's 
%Specify  which  plab.m  PL  coordinates  to  use  and  which  type  of  position 
%(C1,  PR,  smPR,  Clsm,  etc.)  to  calculate 
%Last  Updated  3  Aug  2000 


function  [Clsm,  avg]  =  smooth(PR,  LI,  w,  wvar,  wmin); 


global  file  m  n 
start  =  1 ; 

stop  =  size(PR,l); 
c  =  299792458.0; 
fl  =  1575.42e6; 

Ll_m  =  Ll  *  c/fl; 

Clsm  =  PR; 
once  =  1; 

%Begin  observation  time  loop 
for  t  =  2:m 

C 1  sm(t, :)  =  w*PR(t, :)+( 1  -w)*  (C 1  sm(t- 1 ,  :)+L  1  _m(t ,  :)-L  1  _m(t- 1 , :)) ; 


%speed  of  light  in  m/s 
%L1  frequency  in  MHz 
%  carrier  phase  in  meters 


w  =  w-wvar; 
if  w  <  wmin 
if  once=l 
avg  =  t 
once=2; 
end 

w  =  wmin; 
end 

end 
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